str(dmf2)
class(dmf2$activity)
dmf2 %>% tabyl(activity)

actDmf2 <- filter(dmf2, activity == "TRUE")
actDmf2 %>% tabyl(activity_type)
cookActDmf2 <- filter(actDmf2, activity_type == "cooking")
which(is.na(dmf2$co_pp))
dmf2$co_pp[which(is.na(dmf2$co_pp))] <- 0
dmf2$co_pp <- as.numeric(dmf2$co_pp)
class(dmf2$co_pp)
dmf2 %>% summarise(co_quantiles = quantile(co_pp))
blDmf2 <- dmf2 %>% filter (activity == "FALSE")
blDmf2 %>% summarise(co_quantiles = quantile(co_pp))
actDmf2 %>% summarise(co_quantiles = quantile(co_pp))
cookActDmf2 %>% summarise(co_quantiles = quantile(co_pp))

cookActDmf2 %>% tabyl(Fuel)
cookActDmf2 %>% group_by(Fuel) %>% summarise(medCookPm = median(pm_mean), n=n())
cookActDmf2 %>% group_by(Fuel) %>% summarise(pm_quantiles = quantile(pm_mean))
cookActDmf2 %>% tabyl(Stove)

cookActDmf2 %>% tabyl(Stove, ID)

cookActDmf2 %>% group_by(Stove) %>% summarise(medCookPm = median(pm_mean), n=n())
cookActDmf2 %>% group_by(Stove) %>% summarise(pm_quantiles = quantile(pm_mean))
cookActDmf2 %>% tabyl(Place, ID)

dmf2$Place[dmf2$Place == "full kitchen"] <- "indoors"
dmf2$Place[dmf2$Place == "inside house"] <- "indoors"
actDmf2 <- filter(dmf2, activity == "TRUE")
cookActDmf2 <- filter(actDmf2, activity_type == "cooking")
cookActDmf2 %>% tabyl(Place, ID)
cookActDmf2$placeFctr <- factor(cookActDmf2$Place, order = FALSE)
cookActDmf2 %>% tabyl(placeFctr)
dmf2 %>%
count(ID) %>%
mutate(hr = (n*2)/60)
dmf2 %>%
count(activity) %>%
mutate(hr = (n*2)/60)
dmf2 %>%
count(ID) %>%
mutate(hr = (n*2)/60) %>%
summarise(medTraceH = median(hr))
shapiro.test(cookActDmf2$pm_mean)
actDmf2 %>%
  summarise(quantActpm = quantile(pm_mean))
blDmf2 <- dmf2 %>% filter (activity == "FALSE")
blDmf2 %>% summarise(quantBLpm = quantile(pm_mean))
options(scipen=999)
box_plot_act_dmf2 <- ggplot(dmf2, aes(x = pm_mean, y = activity)) + stat_boxplot() + geom_boxplot() + scale_x_log10() + ggtitle ("Boxplot comparing exposures during combustion activity and baseline") + labs(x = "PM2.5 concentration", y = "Activity/Baseline") + scale_y_discrete(labels = c("Baseline", "Activity"))
box_plot_dnBL_dmf2 <- ggplot(blDmf2, aes(x = pm_mean, y = timeofday)) + geom_boxplot() + scale_x_log10() + ggtitle("Box plot comparing ambient exposures during day and night hours") + labs(x = "PM2.5 concentration", y = "Day/Night")

dayDmf2 <- dmf2 %>% filter (timeofday == "day")
nightDmf2 <- dmf2 %>% filter (timeofday == "night")
dayBlDmf2 <- dayDmf2 %>% filter (activity == "FALSE")
nightBlDmf2 <- nightDmf2 %>% filter(activity == "FALSE")
dayBlDmf2 %>% summarise(dayBaselinepm_quantiles = quantile(pm_mean))
nightBlDmf2 %>% summarise(nightBaselinepm_quantiles = quantile(pm_mean))
cookActDmf2 %>% tabyl(Stove)
cookActDmf2 %>% tabyl(Place)
box_plot_fuel_dmf2 <- ggplot(cookActDmf2, aes(x = pm_mean, y = Fuel)) + geom_boxplot() + scale_x_log10() + labs(y= "Fuel", x = "PM2.5 concentration") + ggtitle("Boxplot comparing cooking-related PM2.5 exposures by cooking fuel")
box_plot_stove_dmf2 <- ggplot(cookActDmf2, aes(x = pm_mean, y = Stove)) + geom_boxplot() + scale_x_log10() + scale_y_discrete(labels = c('three stone fire','charcoal stove','firewood stove')) + labs(y= "Cooking device", x = "PM2.5 concentration") + ggtitle("Boxplot comparing cooking-related PM2.5 exposures by cooking device")
box_plot_place_dmf2 <- ggplot(cookActDmf2, aes(x = pm_mean, y = Place)) + geom_boxplot() + scale_x_log10() + labs(y= "Place of cooking", x = "PM2.5 concentration") + ggtitle("Boxplot comparing cooking-related PM2.5 exposures by place of cooking")
cor.test( ~ log(1+pm_mean) + log(1+co_pp), data = dmf2, method = "spearman", continuity = FALSE, conf.level = 0.95)
library(cowplot)
overall_dens_plot <- ggplot(dmf2, aes(x = log(1+pm_mean), y = log(1+co_pp))) +  geom_density_2d_filled(bins = 9) + scale_y_log10() + xlab("log(1 + PM2.5 concentration)") + ylab("log(1 + CO concentration)") + geom_smooth(method = "lm", colour = "azure2", size = 0.3)
bl_dens_plot <- ggplot(blDmf2, aes(x = log(1+pm_mean), y = log(1+co_pp))) +  geom_density_2d_filled(bins = 9) + scale_y_log10() + xlab("log(1 + PM2.5 concentration)") + ylab("log(1 + CO concentration)") + geom_smooth(method = "lm", colour = "azure2", size = 0.3)
cook_dens_plot <- ggplot(cookActDmf2, aes(x = log(1+pm_mean), y = log(1+co_pp))) +  geom_density_2d_filled(bins = 9) + scale_y_log10() + xlab("log(1 + PM2.5 concentration)") + ylab("log(1 + CO concentration)") + geom_smooth(method = "lm", colour = "azure2", size = 0.3)
cor.test( ~ log(1+pm_mean) + log(1+co_pp), data = blDmf2, method = "spearman", continuity = FALSE, exact = FALSE, conf.level = 0.95)
cor.test( ~ log(1+pm_mean) + log(1+co_pp), data = cookActDmf2, method = "spearman", continuity = FALSE, exact = FALSE, conf.level = 0.95)
noncookDmf2 <- filter(actDmf2, activity_type != "cooking")
cor.test( ~ log(1+pm_mean) + log(1+co_pp), data = noncookDmf2, method = "spearman", continuity = FALSE, exact = FALSE, conf.level = 0.95)
ggplot(noncookDmf2, aes(x = log(1+pm_mean), y = log(1+co_pp))) +  geom_density_2d_filled(bins = 9) + scale_y_log10() + xlab("log(1 + PM2.5 concentration)") + ylab("log(1 + CO concentration)") + geom_smooth(method = "lm", colour = "azure2", size = 0.3)
box_plot_co_fuel_dmf2 <- ggplot(cookActDmf2, aes(x = co_pp, y = Fuel)) + geom_boxplot() + scale_x_log10() + labs(y= "Fuel", x = "CO concentration") + ggtitle("Boxplot comparing cooking-related CO exposures by cooking fuel")
box_plot_co_stove_dmf2 <- ggplot(cookActDmf2, aes(x = co_pp, y = Stove)) + geom_boxplot() + scale_x_log10() + labs(y= "Fuel", x = "CO concentration") + ggtitle("Boxplot comparing cooking-related CO exposures by cooking device")
cookActDmf2 %>% summarise(quantActco = quantile(co_pp))
cookActDmf2 %>% group_by(Fuel) %>% summarise(quantActco = quantile(co_pp))
cookActDmf2 %>% group_by(Stove) %>% summarise(quantActco = quantile(co_pp))
box_plot_co_fuel_dmf2 + geom_vline(xintercept = 3.492, linetype="dashed")
box_plot_co_stove_dmf2 + geom_vline(xintercept = 3.492, linetype="dashed")
box_plot_fuel_dmf2 <- ggplot(cookActDmf2, aes(x = pm_mean, y = Fuel)) + geom_boxplot() + geom_vline(xintercept = 15, linetype="dashed") + scale_x_log10() + labs(y= "Fuel", x = "PM2.5 concentration") + ggtitle("Boxplot comparing cooking-related PM2.5 exposures by cooking fuel")
box_plot_stove_dmf2 <- ggplot(cookActDmf2, aes(x = pm_mean, y = Stove)) + geom_boxplot() + geom_vline(xintercept = 15, linetype="dashed") + scale_x_log10() + scale_y_discrete(labels = c('three stone fire','charcoal stove','firewood stove')) + labs(y= "Cooking device", x = "PM2.5 concentration") + ggtitle("Boxplot comparing cooking-related PM2.5 exposures by cooking device")
box_plot_act_dmf2 <- ggplot(dmf2, aes(x = pm_mean, y = activity)) + stat_boxplot() + geom_boxplot() + geom_vline(xintercept = 15, linetype="dashed") + scale_x_log10() + ggtitle ("Boxplot comparing exposures during combustion activity and baseline") + labs(x = "PM2.5 concentration", y = "Activity/Baseline") + scale_y_discrete(labels = c("Baseline", "Activity"))
box_plot_dnBL_dmf2 <- ggplot(blDmf2, aes(x = pm_mean, y = timeofday)) + geom_boxplot() + geom_vline(xintercept = 15, linetype="dashed") + scale_x_log10() + ggtitle("Box plot comparing ambient exposures during day and night hours") + labs(x = "PM2.5 concentration", y = "Day/Night")
box_plot_co_act_dmf2 <- ggplot(dmf2, aes(x = co_pp, y = activity)) + stat_boxplot() + geom_boxplot() + geom_vline(xintercept = 3.492, linetype="dashed") + scale_x_log10() + ggtitle ("Boxplot comparing CO exposures during combustion activity and baseline") + labs(x = "CO concentration", y = "Activity/Baseline") + scale_y_discrete(labels = c("Baseline", "Activity"))
box_plot_co_dnBL_dmf2 <- ggplot(blDmf2, aes(x = co_pp, y = timeofday)) + geom_boxplot() + geom_vline(xintercept = 3.492, linetype="dashed") + scale_x_log10() + ggtitle("Box plot comparing ambient CO exposures during day and night hours") + labs(x = "CO concentration", y = "Day/Night")
box_plot_co_place_dmf2 <- ggplot(cookActDmf2, aes(x = co_pp, y = Place)) + geom_boxplot() + geom_vline(xintercept = 3.492, linetype="dashed") + scale_x_log10() + labs(y= "Place of cooking", x = "CO concentration") + ggtitle("Boxplot comparing cooking-related CO exposures by place of cooking")
dmf1_fp10 <- dmf1 %>%
  filter(ID == "FP10")
pm_co_trace_fp10 <- ggplot(data = dmf1_fp10, aes(x=newdateTime)) + geom_line( aes(y=pm_mean, colour = "PM2.5 concentration")) + geom_line( aes(y=co_pp * coeff, colour = "CO concentration")) + scale_y_continuous(name = "pm2.5 concentration", sec.axis = sec_axis(~. / coeff, name="CO concentration")) + labs(color = "Legend") + scale_colour_manual(values = colours)
pm_co_trace_fp10

